THE BIERMANN CATASTROPHE IN NUMERICAL MAGNETOHYDRODYNAMICS

As discussed in Drake (2006), a strong shock disturbs the thermal equilibrium between electrons and ions in a plasma. That equilibrium is maintained by electron–ion collisions and operates over timescales ${{\tau }_{ei}}$ that are long compared to the shock-crossing time of a parcel of fluid entering the ion viscous shock in consequence of the large ratio ${{m}_{i}}/{{m}_{e}}$ (Spitzer 1962). As a result, it is essential to describe the fluid in terms of an additional degree of freedom—the electron temperature, T_e—with respect to the usual equilibrium MHD model.

Fortunately, it is unnecessary to model the fluid using new inertial degrees of freedom to describe the electron fluid. A single inertial component for the fluid as a whole gives an adequate description of the fluid structure near the shock and a completely satisfactory description of the fluid in smooth flow regions, where electron–ion collisions restore local thermal equilibrium (Drake 2006; Zel'dovich & Raizer 2002). This makes it easier to adapt existing MHD codes to treat the Biermann effect correctly, since it is much easier to add a scalar degree of freedom for T_e than it would be to deal with two velocity fields, one each for the electron fluid and for the ion fluid.

The large ion-to-electron mass ratio, which, as we recall from the discussion preceding Equation (2), is responsible for the charge separation that produces the Biermann effect, has the further consequence that electron mobility is much higher than ion mobility, and that thermal conductivity due to electron-electron collisions dominates the heat transport in the fluid. At the same time, the forces between electrons and ions as the fluid crosses the ion viscous shock front are effectively dissipation-free, as they are simply electrostatic fields generated by charge separation, and collisional dissipation processes are too slow to act during the shock-crossing timescale.

These observations lead to the conclusion that while the ions undergoing shock compression experience the usual thermodynamically irreversible, entropy-generating process associated with shocks, the electrons are compressed adiabatically by the electrostatic forces exerted upon them by the ions and consequently do not suffer entropy increments due to shock compression. The electron entropy would be a passively advected scalar, then, except for the dissipative effect of electron thermal conduction and for the slow (relative to timescales relevant to the shock) effect of electron–ion collisional equilibration. Electron thermal conduction effectively rules out any sudden change in T_e at the shock, since such a change would produce an enormous restoring heat flux to heal the discontinuity.

The fluid structure that follows from these considerations is of a sudden discontinuous compression of the ions at the shock, accompanied by a smooth increase in T_e, which is continuous throughout the shock (in the Eulerian limit where the shock width tends to zero, T_e acquires a discontinuous derivative in the direction of the shock normal). The electron temperature exhibits a thermal precursor that leads the shock. The size ${{\lambda }_{T}}$ of the precursor region may be calculated by balancing advection against heat diffusion in the frame of the shock and is about

$$\begin{eqnarray}&&{{\lambda }_{T}}\approx \frac{{{\kappa }_{e}}}{\rho {{c}_{v,e}}D},\end{eqnarray} \tag{ 6 }$$

where ${{\kappa }_{e}}$ is the electron thermal conductivity, ${{c}_{v,e}}$ is the electron specific heat at constant volume per unit mass, and D is the shock speed (Shafranov 1957; Zel'dovich & Raizer 2002). The effect of this precursor region on accretion shocks in galaxy clusters has been recently studied in Smith et al. (2013).

The sudden change in the motion of ions entering the shock sheath, together with the more highly mobile motion of the electrons in the vicinity of the shock, results in a zone of charge-separation-driven electric field in the normal direction to the shock (Jaffrin & Probstein 1964; Amendt et al. 2009; Zel'dovich & Raizer 2002). By solving the Navier–Stokes equation across the shock (treating the results as means of kinetic-theory distributions, since a fluid picture is clearly invalid inside the viscous shock sheath), Jaffrin & Probstein (1964) calculated the full electric structure across the shock. We do not require the full structure here, since in the end we wish to adopt an Eulerian picture of the shock, in which the width of the shock is zero. For our purposes, the EMF across the shock is all that is required. It is shown in Amendt et al. (2009) that the EMF $\mathcal{E}$ is given by

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{E} & = & {{k}_{B}}{{T}_{e}}{\rm ln} \left( \frac{{{n}_{e,d}}}{{{n}_{e,u}}} \right) \\ {} & = & {{k}_{B}}{{T}_{e}}{\rm ln} \left( \frac{{{\rho }_{d}}}{{{\rho }_{u}}} \right), \\ \end{array}\end{eqnarray} \tag{ 7 }$$

where the subscripts "u" and "d" denote "upstream" and "downstream," respectively, and where we have assumed complete ionization to obtain the second line.

This electric field is due to charge separation and thus has a common ancestry with the electric field ${{{\boldsymbol{E}} }_{B}}$ in Equation (2). However, it is proper to the shock, rather than to the smooth-flow portion of the fluid. It is essential, for the purpose of physical consistency, to demonstrate that an MHD model implementing the Biermann effect should reproduce the shock-crossing EMF given in Equation (7). That this is possible is demonstrated in the next section.

We wrap up this section by exhibiting the Biermann-effect-laden MHD model implied by these kinetic-theory considerations. The dynamical system describing the model is given in conservation form as follows:

$$\begin{eqnarray}&&\frac{\partial \rho }{\partial t}+\nabla \cdot \left( \rho {\boldsymbol{u}} \right)=0\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\frac{\partial \rho {\boldsymbol{u}} }{\partial t}+\nabla \cdot \left[ \rho {\boldsymbol{uu}} -\frac{1}{4\pi }{\boldsymbol{BB}} +{\boldsymbol{1}} \left( P+\frac{{{{\boldsymbol{B}} }^{2}}}{8\pi } \right) \right]=0\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial }{\partial t}\left[ \rho \left( \frac{1}{2}{{{\boldsymbol{u}} }^{2}}+{{\epsilon }_{T}} \right)+\frac{{{{\boldsymbol{B}} }^{2}}}{8\pi } \right]+\nabla \cdot \left[ \rho {\boldsymbol{u}} \left( \frac{1}{2}{{{\boldsymbol{u}} }^{2}}+{{\epsilon }_{T}}+\frac{P}{\rho } \right)+ \right. \\ \quad \left. \frac{c}{4\pi }\left( -({\boldsymbol{u}} /c)\times {\boldsymbol{B}} +\frac{c\eta }{4\pi }\nabla \times {\boldsymbol{B}} +{{{\boldsymbol{E}} }_{B}} \right)\times {\boldsymbol{B}} \right]=0 \\ \end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial \rho {{s}_{e}}}{\partial t}+\nabla \cdot \left( \rho {\boldsymbol{u}} {{s}_{e}} \right) & = & -T_{e}^{-1}\nabla \cdot \left( -{{\kappa }_{e}}\nabla {{T}_{e}} \right) \\ {} & {} & +\frac{\rho {{c}_{v,e}}}{{{T}_{e}}{{\tau }_{ei}}}({{T}_{i}}-{{T}_{e}}) \\ \end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}+\nabla \cdot \left[ {\boldsymbol{uB}} -{\boldsymbol{Bu}} -\frac{{{c}^{2}}\eta }{4\pi }\nabla {\boldsymbol{B}} \right]+c\nabla \times {{{\boldsymbol{E}} }_{B}}=0.\end{eqnarray} \tag{ 12 }$$

In these equations, $P={{P}_{e}}+{{P}_{i}}$ is the total material pressure, ${{\epsilon }_{T}}$ is the total specific thermal energy per unit mass, and s_e is the specific electron entropy per unit mass.

Equations (8)–(9) are the usual MHD equations of mass and momentum conservation. The energy conservation equation, Equation (10), includes a term to account for resistive dissipation and another to account for the Biermann effect.

Equation (11) expresses the conservation of electron entropy s_e (up to heat conduction and electron–ion heat exchange), which, as we pointed out above, is required by the same approximation ${{m}_{e}}/{{m}_{i}}\to 0$ that gives rise to the Biermann effect in the first place. This equation is introduced to describe the additional degree of freedom required by the two-temperature plasma treatment.

Finally, in the induction equation, Equation (12), the Biermann effect is expressed here as a curl, rather than in the conservation form—a divergence of a flux tensor—that it will assume later.

These dynamical equations are supplemented here by perfect gas equations of state (EOSs) for both the electrons and the ions. We assume total ionization throughout what follows for the sake of simplicity.

By inspection of the Biermann source term in Equation (3), we see that it is proportional to $\nabla {{n}_{e}}\times \nabla {{P}_{e}}$. It follows that field generation by the Biermann effect can only occur if the gradients of n_e and P_e are not aligned. This means that in shocks with planar, cylindrical, or spherical symmetry, the field generation rate is zero. We must therefore treat nonsymmetric shock surfaces in order to analyze nontrivial cases of field generation. Accordingly, in this subsection we set up the basic kinematics of the Biermann effect at general shock surfaces.

To describe the shock surface, we introduce a level function ${\Psi }({\boldsymbol{x}} ,t)$ and use it to define the shock surface ${\Psi }({\boldsymbol{x}} ,t)=0$. Note that the function Ψ is of no dynamical significance, but is rather simply a mathematical convenience for describing the shock surface. We will assume that the shock is moving in the direction of the normal vector ${\boldsymbol{n}} \equiv \nabla {\Psi }/|\nabla {\Psi }|$, so that the region ${\Psi }({\boldsymbol{x}} ,t)\gt 0$ is upstream, whereas the region ${\Psi }({\boldsymbol{x}} ,t)\lt 0$ is downstream. These definitions are illustrated in Figure 1.

Zoom In Zoom Out Reset image size

We denote the local shock speed along ${\boldsymbol{n}}$ by D. By considering the motion of the level surface ${\Psi }({\boldsymbol{x}} ,t)=0,$ it is not difficult to show that

$$\begin{eqnarray}&&D=-\frac{\partial {\Psi }}{\partial t}/\left| \nabla {\Psi } \right|.\end{eqnarray} \tag{ 13 }$$

To describe a moving MHD discontinuity that coincides with the surface ${\Psi }({\boldsymbol{x}} ,t)=0$, we will frequently express a field quantity X by the decomposition $X={{X}_{u}}({\boldsymbol{x}} ,t){\Theta }({\Psi })+{{X}_{d}}({\boldsymbol{x}} ,t){\Theta }(-{\Psi })$, where X_u and X_d are continuous functions, and where ${\Theta }({\Psi })$ is a Heaviside step function. After substituting such expressions into field equations, we will find some terms proportional to the Dirac distribution $\delta ({\Psi })$, resulting from differentiation of the Heaviside functions. We will refer to the collected terms of this form as the "shock part" of the evolution equations. As will be seen, such terms embody R–H jumps, which can thus be efficiently extracted for these nonsymmetric shocks. This trick is also useful for deducing hydrodynamic flux terms, as will also be evident below.

It is convenient to reformulate the Biermann source term using the ideal EOS to replace n_e with T_e, the electron temperature. Assuming an ideal gas EOS, the reformulated electric field is

$$\begin{eqnarray}&&{{{\boldsymbol{E}} }_{B}}\equiv -({{k}_{B}}/e){{T}_{e}}\nabla {\rm ln} {{P}_{e}},\end{eqnarray} \tag{ 14 }$$

so that the source term due to the Biermann effect in the induction equation is

$$\begin{eqnarray}&&{{\left. \frac{\partial {\boldsymbol{B}} }{\partial t} \right|}_{B}}=(c{{k}_{B}}/e)\nabla {{T}_{e}}\times \nabla {\rm ln} {{P}_{e}}.\end{eqnarray} \tag{ 15 }$$

An important reason for this reformulation is that, as we saw in Section 2, in the presence of electron thermal conduction, T_e is continuous at the shock, whereas n_e is not (Shafranov 1957; Zel'dovich & Raizer 2002). As we saw, the continuity of T_e is a consequence of the high mobility of electrons relative to ions and is therefore part and parcel of the same approximation that led to the Biermann source term (Equation (2)) in the first place. It is central to the developments that follow.

The log electron pressure ${\rm ln} {{P}_{e}}$ is discontinuous at the shock and may be represented at time t by

$$\begin{eqnarray}&&{\rm ln} {{P}_{e}}={\rm ln} {{P}_{e,u}}\;{\Theta }\left( {\Psi } \right)+{\rm ln} {{P}_{e,d}}\;{\Theta }\left( -{\Psi } \right),\end{eqnarray} \tag{ 16 }$$

where ${{P}_{e,u}}({\boldsymbol{x}} ,t)$ and ${{P}_{e,d}}({\boldsymbol{x}} ,t)$ are continuous functions.

The discontinuity of P_e leads to a Dirac δ singularity in $\nabla {{P}_{e}}$. Using the distributional relation $d{\Theta }({\Psi })/d{\Psi }=\delta ({\Psi })$, we have

$$\begin{eqnarray}\begin{array}{rcl} \nabla {\rm ln} {{P}_{e}} & = & {\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\delta \left( {\Psi } \right)\nabla {\Psi } \\ {} & {} & +{\Theta }\left( {\Psi } \right)\nabla {\rm ln} {{P}_{e,u}}+{\Theta }\left( -{\Psi } \right)\nabla {\rm ln} {{P}_{e,d}}. \\ \end{array}\end{eqnarray} \tag{ 17 }$$

We may easily interpret the term in Equation (17) that is proportional to $\delta ({\Psi })$ as the gradient proper to the shock, and the remaining terms as the gradient in smooth flow. Since we are interested in the shock behavior, we define

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\left. \nabla {\rm ln} {{P}_{e}} \right|}_{{\rm Shock}}} & \equiv & {\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\delta \left( {\Psi } \right)\nabla {\Psi } \\ {} & = & {\rm ln} \left( \frac{{{\rho }_{u}}}{{{\rho }_{d}}} \right)\delta \left( {\Psi } \right)\nabla {\Psi }, \\ \end{array}\quad \end{eqnarray} \tag{ 18 }$$

where we have used the continuity of T_e and the assumption of complete ionization to replace the pressure ratio across the shock with the corresponding compression ratio ${{\rho }_{u}}/{{\rho }_{d}}$.

The electron temperature T_e is continuous across the shock, but its normal derivative is discontinuous. We may therefore represent T_e by

$$\begin{eqnarray}&&{{T}_{e}}={{T}_{e,u}}{\Theta }\left( {\Psi } \right)+{{T}_{e,d}}{\Theta }\left( -{\Psi } \right),\end{eqnarray} \tag{ 19 }$$

where ${{T}_{e,u}}({\boldsymbol{x}} ,t)$ and ${{T}_{e,d}}({\boldsymbol{x}} ,t)$ are continuous functions, and where continuity at the shock implies

$$\begin{eqnarray}&&{\Psi }({\boldsymbol{x}} ,t)=0\Rightarrow {{T}_{e,u}}({\boldsymbol{x}} ,t)={{T}_{e,d}}({\boldsymbol{x}} ,t).\end{eqnarray} \tag{ 20 }$$

By virtue of the continuity of T_e, the gradient $\nabla {{T}_{e}}$ is not burdened by a Dirac δ, and we obtain

$$\begin{eqnarray}&&\nabla {{T}_{e}}={\Theta }\left( {\Psi } \right)\nabla {{T}_{e,u}}+{\Theta }\left( -{\Psi } \right)\nabla {{T}_{e,d}}.\end{eqnarray} \tag{ 21 }$$

Inserting Equations (18) and (21) into Equation (15), we obtain

$$\begin{eqnarray}\begin{array}{rcl} {{\left. \frac{\partial {\boldsymbol{B}} }{\partial t} \right|}_{B,{\rm Shock}}} & = & (c{{k}_{B}}/e){\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\delta \left( {\Psi }({\boldsymbol{x}} ) \right) \\ {} & {} & \times \left[ \frac{1}{2}\left[ \nabla {{T}_{e,u}}+\nabla {{T}_{e,d}} \right]\times \nabla {\Psi } \right], \\ \end{array}\end{eqnarray} \tag{ 22 }$$

where we have made use of the distributional relation

$$\begin{eqnarray}\begin{array}{rcl} {\Theta }({\Psi })\delta ({\Psi }) & = & {\Theta }({\Psi })\frac{d}{d{\Psi }}{\Theta }({\Psi }) \\ {} & = & \frac{1}{2}\frac{d}{d{\Psi }}{\Theta }{{({\Psi })}^{2}} \\ {} & = & \frac{1}{2}\frac{d}{d{\Psi }}{\Theta }({\Psi }) \\ {} & = & \frac{1}{2}\delta ({\Psi }), \\ \end{array}\end{eqnarray} \tag{ 23 }$$

where, to get the third line, we have used the fact that the Heaviside function squares to itself, since its values are either 0 or 1. This derivation of ${\Theta }\delta =\delta /2$ is admittedly cavalier in its treatment of distributional quantities and is presented for brevity. The result may also be obtained by a limiting procedure, in which the Heaviside function is represented by the limit of a family of continuous functions—this is, after all, the "Eulerian" limit of the description of a shock as the viscosity goes to zero.

Introducing the shock normal vector ${\boldsymbol{n}}$, and using the fact that the tangential derivative of T_e is continuous (so that ${\boldsymbol{n}} \times \nabla {{T}_{e,u}}={\boldsymbol{n}} \times \nabla {{T}_{e,d}}$ at the shock), we obtain

$$\begin{eqnarray}\begin{array}{rcl} {{\left. \frac{\partial {\boldsymbol{B}} }{\partial t} \right|}_{B,{\rm Shock}}} & = & \delta \left( {\Psi } \right)\left| \nabla {\Psi } \right| \\ {} & {} & \times (c{{k}_{B}}/e){\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} . \\ \end{array}\end{eqnarray} \tag{ 24 }$$

Equation (24) is amenable to a clarifying interpretation, which follows from the recognition that if V is a region containing some portion Σ of the shock surface, and the differential element of surface area on Σ is dA, then ${{\int }_{V}}{{d}^{3}}{\boldsymbol{x}} \;\delta ({\Psi }({\boldsymbol{x}} ))|\nabla {\Psi }|={{\int }_{{\Sigma }}}dA$, which is to say that $\delta ({\Psi }({\boldsymbol{x}} ))|\nabla {\Psi }|{{d}^{3}}{\boldsymbol{x}}$ is the differential area element on the shock surface. We may therefore interpret the quantity $(c{{k}_{B}}/e){\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}}$ in Equation (24) as a field generation rate per unit time per unit area of the shock surface. This interpretation already allows us to see that the Biermann source term should give rise to finite, mathematically sensible field generation.

To verify this, we may now obtain the field generation rate due to the passage of the shock by inserting the above expressions into the induction equation, Equation (3), neglecting resistivity. The equation may be written as

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}=\nabla \cdot \left( {\boldsymbol{Bu}} -{\boldsymbol{uB}} \right)+{{\left. \frac{\partial {\boldsymbol{B}} }{\partial t} \right|}_{B}}.\end{eqnarray} \tag{ 25 }$$

In the neighborhood of the shock, we define

$$\begin{eqnarray}&&{\boldsymbol{B}} ={{B}_{n}}{\boldsymbol{n}} +{{{\boldsymbol{B}} }_{Tu}}{\Theta }({\Psi })+{{{\boldsymbol{B}} }_{Td}}{\Theta }(-{\Psi })\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\boldsymbol{u}} ={{{\boldsymbol{u}} }_{u}}{\Theta }({\Psi })+{{{\boldsymbol{u}} }_{d}}{\Theta }(-{\Psi })\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{{\boldsymbol{u}} }_{u}}\equiv {{{\boldsymbol{u}} }_{Tu}}+{{u}_{nu}}{\boldsymbol{n}} \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{{\boldsymbol{u}} }_{d}}\equiv {{{\boldsymbol{u}} }_{Td}}+{{u}_{nd}}{\boldsymbol{n}} ,\end{eqnarray} \tag{ 29 }$$

where the "tangential" components satisfy ${{{\boldsymbol{B}} }_{Tu}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{B}} }_{Td}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{u}} }_{Tu}}\cdot {\boldsymbol{n}} =0$, ${{{\boldsymbol{u}} }_{Td}}\cdot {\boldsymbol{n}} =0$. By virtue of the solenoidal condition $\nabla \cdot {\boldsymbol{B}} =0$, B_n is continuous across the shock.

We transform the advection term $\nabla \cdot \left( {\boldsymbol{Bu}} -{\boldsymbol{uB}} \right)$ in Equation (25) by inserting Equations (26)–(27), grouping the Heaviside functions, and collecting the "shock" terms proportional to the Dirac δ that arises from differentiating the Heaviside functions. Using Equation (13) to replace $\partial {\Psi }/\partial t$, and using Equation (24), we obtain

$$\begin{eqnarray}\begin{array}{rcl} \delta ({\Psi })|\nabla {\Psi }|(-D)({{{\boldsymbol{B}} }_{Tu}}-{{{\boldsymbol{B}} }_{Td}}) & = & \delta ({\Psi })|\nabla {\Psi }| \\ {} & {} & \times \left\{ {{B}_{n}}[{{{\boldsymbol{u}} }_{Tu}}-{{{\boldsymbol{u}} }_{Td}}] \right. \\ \end{array}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\begin{array}{l} -[{{u}_{nu}}{{{\boldsymbol{B}} }_{Tu}}-{{u}_{nd}}{{{\boldsymbol{B}} }_{Td}}] \\ +\left. \frac{c{{k}_{B}}}{e}{\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} \right\}, \\ \end{array}\end{eqnarray} \tag{ 31 }$$

hence the R–H jump condition, generalized by the Biermann flux, is

$$\begin{eqnarray}&&\left[ (D-{{u}_{n}}){{{\boldsymbol{B}} }_{T}} \right]_{d}^{u}+{{B}_{n}}\left[ {{{\boldsymbol{u}} }_{T}} \right]_{d}^{u}+\frac{c{{k}_{B}}}{e}\nabla {{T}_{e}}\times {\boldsymbol{n}} \left[ {\rm ln} {{P}_{e}} \right]_{d}^{u}=0.\end{eqnarray} \tag{ 32 }$$

Here the notation $[\ldots ]_{d}^{u}$ means the difference of the upstream and downstream values at the shock location. The first two terms in Equation (32) constitute the usual jump condition for the induction equation (see, e.g., Chapter 7 of Gurnett & Bhattacharjee 2005). The final term is the contribution to the jump condition from the Biermann effect, which is seen to produce a finite, well-defined discontinuity at the shock. We may obtain a useful result by considering the jump in ${\boldsymbol{B}}$ due to a shock advancing into a quiescent, unmagnetized fluid (B_n = 0, ${{B}_{Tu}}=0$, ${{{\boldsymbol{u}} }_{u}}=0$). The downstream magnetic field is then given by

$$\begin{eqnarray}\begin{array}{rcl} {{{\boldsymbol{B}} }_{Td}} & = & \frac{c{{k}_{B}}}{e}{{(D-{{u}_{nd}})}^{-1}}{\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} \\ {} & = & \frac{c{{k}_{B}}}{eD}\frac{{{\rho }_{d}}}{{{\rho }_{u}}}{\rm ln} \left( \frac{{{\rho }_{u}}}{{{\rho }_{d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} , \\ \end{array}\end{eqnarray} \tag{ 33 }$$

where we have used the condition of mass continuity at the shock, ${{\rho }_{u}}D={{\rho }_{d}}(D-{{u}_{nd}})$, to replace the term $(D-{{u}_{nd}})$, and where we have also used the continuity of T_e to replace the pressure ratio in the log with the density ratio.

We conclude from the above development that the Biermann source term is mathematically well defined even at weak solution discontinuities and yields definite finite predictions to which a properly designed numerical algorithm should be expected to converge.

In the Introduction, we raised the question of whether the Biermann source term, Equation (2), behaves near a shock according to the predictions of kinetic theory, as summarized in Section 2. We now explain the line of thinking behind this question and answer it definitively.

The induction equation, Equation (3), can be cast in conservative form. In this form, the Biermann source term assumes the form of the divergence of a flux tensor, whose components are linear in $\nabla {{P}_{e}}$. It is clear that the construction of these flux components from the gradient of a discontinuous function is in some way associated with the numerical troubles that arise with the Biermann effect near shocks.

The key point here is that "trouble with a flux computed from a derivative" is actually a familiar situation from radiation diffusion, where the radiation flux, which is proportional to the gradient of the radiation temperature, yields unphysical (superluminal) fluxes at regions where the temperature changes sharply (see, e.g., Mihalas & Weibel-Mihalas 1999). Another analogous situation arises with respect to electron thermal conduction, where the thermal flux from $\nabla {{T}_{e}}$ may yield unphysically large transport at shocks, in consequence of the discontinuous change in T_e (Mihalas & Weibel-Mihalas 1999). In both of these cases, a straightforward work-around is furnished by a flux limiter, in which a maximum flux deduced from a kinetic-theory picture of the transport is used to smoothly cut off the flux in regions where gradients get large, while leaving the fluxes in smooth regions unaltered.

This is the reason that the question of the validity of the Biermann flux at shocks is a natural one to ask. If it were the case that the flux is just wrong when $|\nabla {{P}_{e}}|$ gets too large, a reasonable approach would be to use the limiting value for the flux from kinetic theory (Equation 7) and impose that flux in the shock as a cutoff through some kind of flux limiter. If the electric field due to charge separation given by Equation (2) results in an EMF across a shock that exceeds this value, we could cut it off at this value and obtain approximate results analogous to flux-limited approximation-based results from diffusion-limited radiation transport.

Being thus prepared for bad news about the Biermann term from kinetic theory, we instead are met by a surprise upon examination of the MHD flux. The Biermann electric field at a shock is, according to Equations (14) and (18),

$$\begin{eqnarray}&&{{{\boldsymbol{E}} }_{B}}=({{k}_{B}}/e){{T}_{e}}{\rm ln} \left( \frac{{{\rho }_{d}}({\boldsymbol{x}} )}{{{\rho }_{u}}({\boldsymbol{x}} )} \right)\delta \left( {\Psi }({\boldsymbol{x}} ) \right)\nabla {\Psi }.\end{eqnarray} \tag{ 34 }$$

Integrating this over a vanishingly short shock-crossing path l along the normal to the shock, we obtain the electromagnetic work done by an electron in the presence of the Biermann electric field,

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{E}}_{B}} & = & e{{\int }_{l}}d{\boldsymbol{x}} \cdot {{{\boldsymbol{E}} }_{B}} \\ {} & = & {{k}_{B}}{{T}_{e}}{\rm ln} \left( \frac{{{\rho }_{d}}}{{{\rho }_{u}}} \right), \\ \end{array}\end{eqnarray} \tag{ 35 }$$

which is the same as the value inferred from kinetic theory, Equation (7).

We conclude from this argument that there is no analogy to flux-limited diffusion in the behavior of the Biermann source term at shocks. The "plain vanilla" source gives the correct EMF across the shock and needs no flux limiter to cut it off there. It should be perfectly possible to construct a physically valid algorithm to represent the Biermann MHD source term unmoderated by limiters, one that faithfully reproduces the predictions of kinetic theory at a shock.

As discussed in the Introduction, MHD simulations incorporating the Biermann battery effect have invariably produced catastrophic nonconvergent behavior as soon as discontinuities develop. There are two (related) kinds of catastrophes that have manifested themselves: simulations with spherical or planar symmetry, which develop a nonzero field at shocks with a field intensity that grows with increasing resolution, despite the fact that the charge-separation electric field is irrotational and the magnetic field generation rate should therefore be zero (Fatenejad et al. 2013); and simulations of HEDP laser experiments containing spatially inhomogeneous shocks, similar to the (highly simplified) simulations presented below in Section 6.2.2, in which the field generated at shocks, which is not expected to be zero, fails to converge to a finite value with increasing resolution—in order to obtain reasonable results at all in such simulations, the cells participating in shocks, contact, and material discontinuities must be located at each time step, and the generation rate in those cells artificially set to zero (Fatenejad et al. 2013). Such simulation results are clearly not correct, since they do not treat magnetogeneration due to the ion viscous shock, but at least they do converge with resolution.

Clearly something goes very wrong with the numerics when the flow develops discontinuous behavior, whereas the behavior of the Biermann term in smooth flow is finite, convergent, and stable. To investigate the origin of these catastrophes, we have addressed the following question: suppose the electric field ${{{\boldsymbol{E}} }_{B}}$ of Equation (15) is irrotational, in the sense that

$$\begin{eqnarray}&&\nabla \times {{{\boldsymbol{E}} }_{B}}=\frac{{{k}_{B}}}{e}\nabla {{T}_{e}}\times \nabla {\rm ln} {{P}_{e}}=0.\end{eqnarray} \tag{ 36 }$$

To what extent is this irrotational character preserved after T_e and P_e are separately discretized in a volume-based scheme?

We consider the discretization of two scalar functions $f({\boldsymbol{x}} )$ and $g({\boldsymbol{x}} )$ that are assumed to have gradients that are everywhere collinear, so that there exists some function $\alpha ({\boldsymbol{x}} )$ such that $\nabla g=\alpha \nabla f$. A pair of such functions satisfy $\nabla g\times \nabla f=0$, so that the vector field $g\nabla f$ is irrotational. The Taylor series of such a pair of functions are related to each other by the collinearity constraint. We Taylor-expand both functions to third order, average the resulting expansion over a grid of control volumes, and compute the "discrete curl"—that is, the circulation integral about a cell face—of the discretized "$g\nabla f$" vector field. We omit the details of this extremely lengthy calculation for the sake of brevity and merely present the conclusions that may be drawn therefrom.

The result of the calculation is that the leading-order behavior of the discrete curl of such a discretized vector field is $\mathcal{O}({{{\Delta }}^{2}})$, where Δ is the grid spacing. The coefficient multiplying this term is homogeneous of the order of 3 in the derivatives of f and g. What this means is that in regions of space where the two functions are smooth (C² or smoother), the coefficient has a finite approximation in the discretization, and the discrete curl tends to zero as ${{{\Delta }}^{2}}$ with increasingly fine resolution, yielding the correct curl of an irrotational vector field in the limit.

On the other hand, in the presence of a discontinuity, the coefficient of the discrete curl has no finite approximation in the discretization, but rather diverges as ${{{\Delta }}^{-3}}$ when ${\Delta }\to 0$, reflecting the meaninglessness of the Taylor approximation in proximity to a discontinuity. As a consequence, the discrete curl as a whole diverges as ${{{\Delta }}^{-1}}$ with vanishing Δ in the neighborhood of a discontinuity.

Setting $f={\rm ln} {{P}_{e}}$ and $g={{T}_{e}}$, we recognize immediately the source of the discrete pathology described above. The discretization of P_e and T_e produces a spurious solenoidal component of ${{{\boldsymbol{E}} }_{B}}$ above and beyond any real, physically correct solenoidal component. This solenoidal anomaly is small and converges to zero with increasing resolution in smooth flow. Near a discontinuity, however, the anomaly is not convergent, but rather diverges as ${{{\Delta }}^{-1}}$.

This is the explanation of the numerical Biermann catastrophe. It is, fundamentally, a completely predictable failure of naive, stencil-based approximations to the Biermann source term (Equation (2)), which are not meaningful when Taylor-series approximations to the fluid variables fail in the presence of a discontinuity. The failure is irreducible and does not depend on whether the Biermann flux is added to other MHD fluxes or computed separately as a split term.

The above diagnosis of the origin of the numerical Biermann catastrophe does not directly suggest a cure for the problem. That a cure should exist, however, is strongly suggested by the fact that the analytic treatment of the Biermann effect at shocks in Section 3.1 leads to finite and physically sensible field generation rates. We therefore now partly retrace that development with a view to casting those results in a form suitable for formulation as a numerical algorithm to be incorporated in a Godunov-type volume-based MHD scheme.

In order to clarify this proposed strategy, let us briefly review how volume-based hydrodynamics/MHD schemes work (for details, see, e.g., LeVeque 2002). Nonlinear systems of hyperbolic partial differential equations (PDEs) such as MHD are prone to develop discontinuities such as shocks and contacts. Direct finite-difference discretization approaches break down when this occurs. This breakdown is circumvented by appealing to the conservative nature of the equations and replacing the differential equations with their corresponding conservation-integral forms, applied to a regular mesh of control volumes. Field quantities are interpreted as volume averages over control volumes and are updated by time-centered (i.e., half-time-step forward) fluxes across control-volume faces. These fluxes are finite even when solution discontinuities cause the corresponding differential PDE source terms to misbehave.

The crucial procedure for carrying out such schemes successfully is the computation of the time-centered fluxes by means of the solution of Riemann shock tube problems at the cell interfaces. A Riemann solver takes as its input piecewise-constant initial data with a discontinuity at the interface and calculates the resulting propagating state, comprising piecewise-constant regions separated by discontinuity waves and rarefactions belonging to different families of characteristics. This propagating state is used to infer the MHD variables at the interface a half time step following the initial state. MHD fluxes computed from these variables are then used to update the cell-averaged field quantities (LeVeque 2002; Toro 2009).

Riemann solvers relate adjoining MHD states separated by a discontinuity traveling with wave speed D by means of the R–H relations,

$$\begin{eqnarray}&&D\times ({{{\bf \Phi }}_{{\rm L}}}-{{{\bf \Phi }}_{{\rm R}}})={\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm L}}})-{\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm R}}}).\end{eqnarray} \tag{ 37 }$$

Here ${\bf \Phi }$ is the state vector of conserved field quantities, ${{{\bf \Phi }}^{T}}\equiv [\rho ,\rho {{{\boldsymbol{u}} }^{T}},\rho \mathcal{E},{{{\boldsymbol{B}} }^{T}}]$, that is, the mass, momentum, and energy densities and the magnetic field. ${\boldsymbol{F}}$ denotes the vector of fluxes corresponding to ${\bf \Phi }$. The subscripts "L" and "R" denote "left" and "right" states, respectively, and D is positive for a wave traveling from left to right. (Toro 2009).

We can now state a minimal requirement in order to assimilate the Biermann flux to the other MHD fluxes used to update the fluid state: We need to establish the terms added by the Biermann effect to the flux vector ${\boldsymbol{F}}$. Once these additive fluxes are determined, they can be added to the ideal MHD fluxes.

The additional flux of ${\boldsymbol{B}}$ due to the Biermann effect is obtainable by solving for the shock part of the "Biermann-only" induction equation, Equation (14). That is to say, we assume the distributional forms of Equations (16), (19), and (26), plug them into Equation (15), and keep only the terms proportional to the Dirac δ. The result may be obtained by inspection of Equation (32), setting the velocities ${\boldsymbol{u}}$ to zero:

$$\begin{eqnarray}&&D({{{\boldsymbol{B}} }_{Tu}}-{{{\boldsymbol{B}} }_{Td}})=-\frac{c{{k}_{B}}}{e}{\rm ln} \left( \frac{{{P}_{e,u}}}{{{P}_{e,d}}} \right)\nabla {{T}_{e}}\times {\boldsymbol{n}} ,\end{eqnarray} \tag{ 38 }$$

Comparing Equations (38) and (37) and using the continuity of $\nabla {{T}_{e}}\times {\boldsymbol{n}}$, we obtain for the Biermann effect flux of ${\boldsymbol{B}}$ along ${\boldsymbol{n}}$

$$\begin{eqnarray}&&{{{\boldsymbol{f}} }_{{\boldsymbol{B}} }}({\boldsymbol{n}} )\equiv -\frac{c{{k}_{B}}}{e}{\rm ln} ({{P}_{e}})\nabla {{T}_{e}}\times {\boldsymbol{n}} .\end{eqnarray} \tag{ 39 }$$

This expression is manifestly well defined at discontinuities. As in Equation (32), from which Equation (39) is derived, the only residue of the previously toxic gradient $\nabla {{P}_{e}}$ is contained in the benign direction vector ${\boldsymbol{n}}$, which is effectively $-\nabla {{P}_{e}}/|\nabla {{P}_{e}}|$ at the shock. The effect of the cross product ${\boldsymbol{n}} \times \nabla {{T}_{e}}$ is to eliminate the normal component of $\nabla {{T}_{e}}$, while leaving only the tangential components. It is clear from the form of Equation (39) that only the tangential components of ${\boldsymbol{B}}$ are subject to change according to the R–H condition, as expected.

It is convenient for the purposes of algorithmic implementation to re-express Equation (39) in tensor form, as an antisymmetric tensor ${{{\boldsymbol{G}} }^{(B)}}$ whose components $G_{ik}^{(B)}$ express the flux in the coordinate direction k of the magnetic field component B_i:

$$\begin{eqnarray}&&G_{ik}^{(B)}\equiv -\frac{c{{k}_{B}}}{e}{\rm ln} ({{P}_{e}})\mathop{\sum }\limits_{j=1}^{3}{{\epsilon }_{ijk}}\frac{\partial {{T}_{e}}}{\partial {{x}_{j}}},\end{eqnarray} \tag{ 40 }$$

where ${{\epsilon }_{ijk}}$ is the totally antisymmetric Levi–Civita tensor. $G_{ik}^{(B)}$ is, of course, the spatial part of the dual Maxwell tensor due to the Biermann effect. It is worth pointing out that this expression for the flux tensor differs from the "naive" flux

$$\begin{eqnarray}&&G_{ik}^{(B,{\rm Naive})}\equiv +\frac{c{{k}_{B}}}{e}{{T}_{e}}\mathop{\sum }\limits_{j=1}^{3}{{\epsilon }_{ijk}}\frac{\partial {\rm ln} ({{P}_{e}})}{\partial {{x}_{j}}},\end{eqnarray} \tag{ 41 }$$

obtainable directly from the Biermann electric field (Equation (14)) by a pure gradient, which has no effect on the induction equation. It follows that Equation (39) yields a general flux of ${\boldsymbol{B}}$ in the direction ${\boldsymbol{n}}$ that can be used anywhere in the fluid, not only at discontinuities.

The flux in Equation (39) is not the only correction that must be applied to ${\boldsymbol{F}}$. There is also a correction to be applied to the energy part of the flux vector, by virtue of the Poynting flux $(c/4\pi ){{{\boldsymbol{E}} }_{B}}\times {\boldsymbol{B}}$ that arises in connection with the charge-separation electric field ${{{\boldsymbol{E}} }_{B}}$. This term is a bit puzzling at first: since ${{{\boldsymbol{E}} }_{B}}$ is normal to the discontinuity, the Poynting flux is tangential to the discontinuity, and it is not immediately obvious how such a flux is to be pressed into service in an R–H relation.

Nonetheless, the required flux may be inferred in a manner analogous to the calculation by which we derived ${{{\boldsymbol{F}} }_{{\boldsymbol{B}} }}$, above: We solve the "Biermann-only" energy equation

$$\begin{eqnarray}&&\frac{\partial (\rho \mathcal{E})}{\partial t}+\nabla \cdot \left( \frac{c}{4\pi }{{{\boldsymbol{E}} }_{B}}\times {\boldsymbol{B}} \right)=0,\end{eqnarray} \tag{ 42 }$$

inserting the distributional forms of Equations (16), (19), (26), as well as

$$\begin{eqnarray}&&(\rho \mathcal{E})({\boldsymbol{x}} ,t)={{(\rho \mathcal{E})}_{u}}{\Theta }\left( {\Psi } \right)+{{(\rho \mathcal{E})}_{d}}{\Theta }\left( -{\Psi } \right).\end{eqnarray} \tag{ 43 }$$

After collecting terms proportional to the Dirac δ, we obtain

$$\begin{eqnarray}\begin{array}{rcl} D\left[ \rho \mathcal{E} \right]_{d}^{u} & = & \frac{c{{k}_{B}}}{4\pi e}\left[ \left( ({{T}_{e}}\nabla {\rm ln} {{P}_{e}})\times {\boldsymbol{n}} \right)\cdot {\boldsymbol{B}} \right. \\ {} & {} & -\left( \nabla ({{T}_{e}}{\rm ln} {{P}_{e}})\times {\boldsymbol{n}} \right)\cdot \lt {\boldsymbol{B}} \gt \\ {} & {} & +\left. {{T}_{e}}{\rm ln} {{P}_{e}}{\boldsymbol{n}} \cdot \nabla \times \lt {\boldsymbol{B}} \gt \right]_{d}^{u}, \\ \end{array}\end{eqnarray} \tag{ 44 }$$

where in Equation (44) we distinguish between ${\boldsymbol{B}}$—the magnetic field to the left or right of the interface—and $\lt {\boldsymbol{B}} \gt \equiv \frac{1}{2}({{{\boldsymbol{B}} }_{u}}+{{{\boldsymbol{B}} }_{d}})$, the average of the upstream and downstream fields at the interface. Comparing Equation (44) and Equation (37), we obtain for the energy flux along ${\boldsymbol{n}}$ due to the Biermann effect

$$\begin{eqnarray}\begin{array}{rcl} {{f}_{\rho \mathcal{E}}}({\boldsymbol{n}} ) & \equiv & \frac{c{{k}_{B}}}{4\pi e}\left[ \left( ({{T}_{e}}\nabla {\rm ln} {{P}_{e}})\times {\boldsymbol{n}} \right)\cdot {\boldsymbol{B}} \right. \\ {} & {} & -\left( \nabla ({{T}_{e}}{\rm ln} {{P}_{e}})\times {\boldsymbol{n}} \right)\cdot \lt {\boldsymbol{B}} \gt \\ {} & {} & +\left. {{T}_{e}}{\rm ln} {{P}_{e}}\;{\boldsymbol{n}} \cdot \nabla \times \lt {\boldsymbol{B}} \gt \right]. \\ \end{array}\end{eqnarray} \tag{ 45 }$$

Here we see a potential problem: the third term in Equation (45) requires an average of the upstream and downstream currents $\nabla \times {\boldsymbol{B}}$. This is a difficulty if a true flux is to be extracted and pressed into service in a numerical scheme, since possible discontinuities in ${\boldsymbol{B}}$ make such a term ambiguous. The resolution of the difficulty is that, in general, there is always some finite resistivity in real plasmas. As we will see below in Section 4.1, the presence of finite resistivity removes the discontinuity in ${\boldsymbol{B}} ,$ allowing us to replace $\lt {\boldsymbol{B}} \gt$ with ${\boldsymbol{B}}$. We then obtain for the flux of energy in the direction ${\boldsymbol{n}}$

$$\begin{eqnarray}&&{{f}_{\rho \mathcal{E}}}({\boldsymbol{n}} )=\frac{c{{k}_{B}}}{4\pi e}{\rm ln} {{P}_{e}}\left[ \nabla \times \left( {{T}_{e}}{\boldsymbol{B}} \right) \right]\cdot {\boldsymbol{n}} .\end{eqnarray} \tag{ 46 }$$

To summarize, at a cell interface with a normal vector ${\boldsymbol{n}}$, the hydrodynamic fluxes are to be adapted to the Biermann effect by adding to the implemented flux vector ${\boldsymbol{F}}$ an additional flux vector ${{{\boldsymbol{F}} }^{(B)}}({\boldsymbol{n}} )$, given by

$$\begin{eqnarray}&&{{{\boldsymbol{F}} }^{(B)}}=\left( \begin{array}{ccccccccccccccc} F_{\rho }^{(B)} \\ {\boldsymbol{F}} _{\rho {\boldsymbol{u}} }^{(B)} \\ F_{\rho \mathcal{E}}^{(B)} \\ {\boldsymbol{F}} _{{\boldsymbol{B}} }^{(B)} \\ \end{array} \right)\equiv \left( \begin{array}{ccccccccccccccc} 0 \\ {\boldsymbol{0}} \\ {{f}_{\rho \mathcal{E}}}({\boldsymbol{n}} ) \\ {{{\boldsymbol{f}} }_{{\boldsymbol{B}} }}({\boldsymbol{n}} ) \\ \end{array} \right).\end{eqnarray} \tag{ 47 }$$

We expect a discretization of the Biermann effect based on these flux expressions to be convergent so long as the mesh resolves the scales on which T_e and ${\boldsymbol{B}}$ are continuous. The scale at which T_e is continuous is the scale length ${{\lambda }_{T}}$ of the electron thermal conduction precursor region (Shafranov 1957; Jaffrin & Probstein 1964; Zel'dovich & Raizer 2002), discussed in Section 2.2. This length scale is characteristic of the variation of temperature perpendicular to the shock; however, it is also the length scale over which heat diffuses transversely during the time required for the shock to travel a distance ${{\lambda }_{T}}$ (and hence traverse its own thermal precursor). It is therefore the scale to be resolved in order for the discretization of Equation (39) to converge. At coarser resolutions, T_e appears as discontinuous at shocks as all the other fluid variables, and the discretization discussed here will not yield converged results for ${\boldsymbol{B}}$. Only as the thermal precursor zone is resolved can convergence be expected.

In order for the discretization to produce convergent results for energy, it is necessary that the simulation resolve the resistive length scale discussed in Section 4.1. Note, however, that for many physical situations of interest, the magnetic field intensities generated by the Biermann effect are so tiny that the correction due to the Biermann effect to the energy flux may justifiably be neglected—it is dwarfed by the fluxes of thermal and kinetic energy, and possibly also by the advected flux of existing magnetic field energy. Under such circumstances, it is an acceptable approximation to neglect the Biermann energy flux, if a non-resistive calculation is of interest, or if the resistive scale cannot be resolved.

When the resistivity η in the induction equation, Equation (3), is finite, a new effect appears in the vicinity of a discontinuity: a resistive magnetic precursor travels ahead of the discontinuity, as the impulsively generated field due to the Biermann effect at the shock diffuses out. This effect is analogous to the well-known electron thermal precursor that precedes a plasma shock, which was discussed in Section 2.2. The structure of the thermal precursor can be estimated by balancing diffusion and advection of thermal energy in the frame of the shock, in the presence of impulsive shock heating (Zel'dovich & Raizer 2002). Similarly, as we will presently see, the structure of the magnetic precursor is a consequence of the balance between magnetic field advection and diffusion in the frame of the shock, in the presence of impulsive magnetogeneration at the shock.

We consider a small portion of the shock surface, which we will treat as approximately planar, and we will work in the rest frame of the shock. We assume an approximately steady state in the shock frame, so that we may set $\partial /\partial t=0$ in the field equations. The induction equation then becomes

$$\begin{eqnarray}&&\begin{array}{l} \nabla \cdot \left( {\boldsymbol{B}} {{{\boldsymbol{u}} }^{\prime }}-{{{\boldsymbol{u}} }^{\prime }}{\boldsymbol{B}} \right)+\nabla \cdot \left( \frac{{{c}^{2}}\eta }{4\pi }\nabla {\boldsymbol{B}} \right) \\ \quad +\frac{c{{k}_{B}}}{e}\nabla {{T}_{e}}\times \nabla {\rm ln} {{P}_{e}}={\boldsymbol{0}} , \\ \end{array}\end{eqnarray} \tag{ 48 }$$

where ${{{\boldsymbol{u}} }^{\prime }}$ is the fluid velocity in the rest frame of the shock.

We will keep only the shock part of the Biermann flux in Equation (48), that is, the right-hand side of Equation (24). In effect, this amounts to assuming that the magnetic field generation from the Biermann effect is impulsive at the shock, so that the smooth-flow Biermann effect is much smaller than the effect at the shock. This approximation is justified if the size ${{\lambda }_{T}}$ of the electron thermal precursor is much smaller than the size of ${{\lambda }_{B}}$ of the magnetic precursor, so that thermal gradients are unimportant except at the discontinuity. The plausibility of this condition will be verified in Section 4.3.

We may adopt the local simplification ${\Psi }({\boldsymbol{x}} ,t)={\boldsymbol{n}} \cdot {\boldsymbol{x}}$, $|\nabla {\Psi }|=|{\boldsymbol{n}} |=1$ for the discontinuity-tracing level function. We further adopt coordinates such that x₁ is along ${\boldsymbol{n}}$, so that ${\boldsymbol{n}} ={{{\boldsymbol{e}} }_{1}}$, and such that x₂ is along $\nabla {{T}_{e}}\times {\boldsymbol{n}}$. Away from the discontinuity, we ignore all derivatives except for $\partial /\partial {{x}_{1}}$.

The presence of the resistive term changes the discontinuous structure of the solution. This term is the divergence of the resistive flux of ${\boldsymbol{B}}$. That flux behaves analogously to the Fick's law heat flux $-\kappa \nabla T$, in that it opposes gradients in ${\boldsymbol{B}}$. A discontinuous ${\boldsymbol{B}}$ is disallowed, because it would result in an infinite restoring resistive flux. ${\boldsymbol{B}}$ is therefore now continuous. By the R–H condition for transverse momentum (Gurnett & Bhattacharjee 2005, Chapter 7), it follows that ${\boldsymbol{u}} _{T}^{\prime }$ is also continuous. Only $u_{n}^{\prime }$ and P_e have solution discontinuities in this case. Assuming that the upstream fluid is at rest in the lab frame, we therefore set ${{{\boldsymbol{u}} }^{\prime }}=u_{n}^{\prime }{{{\boldsymbol{e}} }_{1}}$, with

$$\begin{eqnarray}&&u_{n}^{\prime }=(-D){\Theta }({{x}_{1}})+({{u}_{nd}}-D){\Theta }(-{{x}_{1}}).\end{eqnarray} \tag{ 49 }$$

By the coordinate choice, assuming that the far upstream fluid is unmagnetized, and by virtue of Equation (48), we may set ${\boldsymbol{B}} =B{{{\boldsymbol{e}} }_{2}}$. Putting this all together, we obtain the following structure equation for the magnetic precursor:

$$\begin{eqnarray}&&\begin{array}{l} \frac{d}{d{{x}_{1}}}\left[ \frac{{{c}^{2}}\eta }{4\pi }\frac{dB}{d{{x}_{1}}}-u_{n}^{\prime }B \right]+\delta ({{x}_{1}})\frac{c{{k}_{B}}}{e}|\nabla {{T}_{e}}\times {\boldsymbol{n}} |\left[ {\rm ln} {{P}_{e}} \right]_{d}^{u} \\ \quad =0. \\ \end{array}\end{eqnarray} \tag{ 50 }$$

The meaning of this equation is that the precursor structure is determined by the balance of magnetic diffusion and magnetic advection, in the presence of impulsive magnetic field generation due to the Biermann effect.

In the upstream (${{x}_{1}}\gt 0$) region, we may neglect gradients in η (which depends on T_e). Equation (50) becomes

$$\begin{eqnarray}&&\frac{{{c}^{2}}\eta }{4\pi }\frac{{{d}^{2}}B}{dx_{1}^{2}}+D\frac{dB}{d{{x}_{1}}}=0,\end{eqnarray} \tag{ 51 }$$

the solution of which, given ${{{\rm lim} }_{{{x}_{1}}\to +\infty }}B({{x}_{1}})=0,$ is

$$\begin{eqnarray}\begin{array}{rcl} {{B}_{u}}({{x}_{1}}) & = & {{B}_{0}}{\rm exp} \left( -\frac{4\pi D}{{{c}^{2}}\eta }{{x}_{1}} \right) \\ {} & = & {{B}_{0}}{\rm exp} \left( -\frac{{{x}_{1}}}{{{\lambda }_{B}}} \right),\quad {{x}_{1}}\gt 0. \\ \end{array}\end{eqnarray} \tag{ 52 }$$

Here B₀ is an integration constant. We see that the precursor has an exponential shape and a characteristic length ${{\lambda }_{B}}$ given by

$$\begin{eqnarray}&&{{\lambda }_{B}}=\frac{{{c}^{2}}\eta }{4\pi D}.\end{eqnarray} \tag{ 53 }$$

Here $\eta =\eta ({{T}_{e,+\infty }})$ is the value of the resistivity far upstream of the discontinuity.

We may obtain a relation for B₀ from the jump condition at the discontinuity, by integrating Equation (50) across a vanishingly small shock crossing path. The result is

$$\begin{eqnarray}&&\begin{array}{l} \frac{{{c}^{2}}\eta ({{T}_{e}})}{4\pi }\left[ \frac{dB}{d{{x}_{1}}} \right]_{d}^{u}-\left[ {{u}_{n}} \right]_{d}^{u}{{B}_{0}} \\ \quad +\frac{c{{k}_{B}}}{e}|\nabla {{T}_{e}}\times {\boldsymbol{n}} |\left[ {\rm ln} {{P}_{e}} \right]_{d}^{u}=0, \\ \end{array}\end{eqnarray} \tag{ 54 }$$

where now $\eta ({{T}_{e}})$ is evaluated at the shock. This condition, together with the upstream solution and boundary condition, determines B₀.

The ideal MHD jump conditions may be recovered in the limit of spatially constant η with $\eta \to 0$. In that case, the downstream solution satisfying finite-B boundary conditions at ${{x}_{1}}\to -\infty$ is ${{B}_{d}}({{x}_{1}})={{B}_{0}}$, and the first term in Equation (54) is just $-D{{B}_{0}}$. Equation (33) follows immediately.

One further consequence of the presence of the thermal precursor in T_e described in Section 2.2 is that, in general, there will be some small amount of magnetic field generated by the smooth-flow Biermann effect ahead of the shock, owing to plasma motions in the preheated region, whose size is ${{\lambda }_{T}}$, given by Equation (6). In general, the field intensity due to this effect can be expected to be small compared to the field intensity due to the resistive magnetic precursor, since the gradients in the thermal magnetic precursor are small compared to those available near the shock itself. In the constant-conductivity simulations that we discuss in Section 6.2.4, we will see that the field intensity is in fact quite small.

Nonetheless, this smallness does not necessarily preclude the possibility that the thermal magnetic precursor might be experimentally observable. Thermal conductivity upstream of a plasma shock is, in general, not constant, but rather depends strongly on electron temperature—${{\kappa }_{e}}\sim T_{e}^{5/2}$ (Spitzer 1962). The actual structure of the thermal precursor is not the gentle exponential decay that one expects for a constant ${{\kappa }_{e}}$, but rather exhibits the steep gradient near its outer terminus shown in Figure 7.20 of Zel'dovich & Raizer (2002). It is possible that this large gradient may be responsible for an enhancement of the Biermann effect at the terminus of the thermal conduction precursor zone. This interesting possibility is beyond the scope of this paper.

One consequence of the presence of a thermal magnetic precursor is that even in a nonresistive approximation, the notion of an unmagnetized fluid upstream of the shock, on which, for example, Equation (33) is based, is not strictly correct. It is a valid approximation only when the Biermann generation rate in the thermal precursor is small compared to that of the shock.

It is useful to establish some of the relevant length scales under conditions of interest. Below, we establish these scales for conditions prevailing in HEDP experiments. In Section 5 we will establish the scales characteristic of galaxy cluster formation.

Using the expression for η from Huba (2007) in the expression for the resistive magnetic precursor length ${{\lambda }_{B}}$ given in Equation (53), we have the following expression for ${{\lambda }_{B}}$ in a fully ionized plasma:

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{B}} & = & (4/3){{(2\pi )}^{-1/2}}{{c}^{2}}{{e}^{2}}m_{e}^{1/2}Z{\rm ln} {\Lambda }{{D}^{-1}}{{({{k}_{b}}{{T}_{e}})}^{-3/2}} \\ {} & = & 16.4\;{\rm cm}\times Z\times \frac{{\rm ln} {\Lambda }}{10} \\ {} & {} & \times {{\left( \frac{D}{{{10}^{6}}\;{\rm cm}{{{\rm s}}^{-1}}} \right)}^{-1}}\times {{\left( \frac{{{k}_{B}}{{T}_{e}}}{1{\rm eV}} \right)}^{-3/2}}, \\ \end{array}\end{eqnarray} \tag{ 55 }$$

where Λ is the usual Coulomb logarithm. Similarly, the expression for the electron thermal precursor length scale ${{\lambda }_{T}}$ is (Zel'dovich & Raizer 2002, Chapter 7, Section 12)

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{T}} & = & \frac{2}{5}\frac{{{\kappa }_{e}}({{T}_{e,0}})}{\rho {{c}_{v,e}}D} \\ {} & = & \frac{4}{15}\frac{\xi {{({{k}_{B}}{{T}_{e}})}^{5/2}}}{{{e}^{4}}m_{e}^{1/2}ZA{\rm ln} {\Lambda }{{n}_{i}}D} \\ {} & = & 5.39\times {{10}^{-3}}{\rm cm}\times \xi \times {{Z}^{-1}}{{A}^{-1}}\times {{\left( \frac{{\rm ln} {\Lambda }}{10} \right)}^{-1}} \\ {} & {} & \times {{\left( \frac{D}{{{10}^{6}}\;{\rm cm}{{{\rm s}}^{-1}}} \right)}^{-1}}{{\left( \frac{{{k}_{B}}{{T}_{e}}}{1{\rm eV}} \right)}^{5/2}}{{\left( \frac{{{n}_{i}}}{{{10}^{16}}{\rm c}{{{\rm m}}^{-3}}} \right)}^{-1}}, \\ \end{array}\end{eqnarray} \tag{ 56 }$$

where ${{T}_{e,0}}$ is the electron temperature at the shock, and where ξ is a number in the range 1–2.

The numbers in Equations (55) and (56) have been scaled to conditions that are routinely obtainable in large laser facilities such as Omega and NIF (see, e.g., Gregori et al. 2012). In these conditions, it is clear that the assumption ${{\lambda }_{T}}\ll {{\lambda }_{B}}$, required for the validity of the derivation of ${{\lambda }_{B}}$, is easily satisfied.

It is also immediately apparent that the magnetic precursor length scale is a macroscopic scale that can, in principle, be well tailored to the physical dimensions of the interaction chambers of such facilities. A carefully designed experiment, which launches a shock into a relatively cold upstream plasma (note the dependence of ${{\lambda }_{B}}$ on T_e) should in principle be able to detect the resistive precursor due to the Biermann effect.

Note that as the value of T_e upstream rises, ${{\lambda }_{B}}$ decreases as $T_{e}^{-3/2}$, while ${{\lambda }_{T}}$ increases as $T_{e}^{5/2}$. We therefore only need to increase the upstream temperature to about 7 eV for ${{\lambda }_{T}}$ to become comparable to ${{\lambda }_{B}}$, and at warmer temperatures than this ${{\lambda }_{T}}$ becomes dominant. It is therefore plausible that a physical situation could be created in which the thermal magnetic precursor might be observable without interference from the resistive magnetic precursor.

As discussed in the Introduction, the Biermann effect has been invoked as the source of seed magnetic fields that can be amplified by the turbulent dynamo mechanism (Kulsrud et al. 1997).

Given the fact that in these studies the gradients near shocks that were used to calculate the Biermann effect field generation rates are artifacts of the hydrodynamic advance schemes, it is possible—even likely—that the calculated field strengths are incorrect. We now attempt to estimate the typical strength of magnetic fields that are to be expected in early galaxy cluster formation, based on a corrected treatment of the Biermann effect. We do not perform new simulations, but rather attempt to estimate the typical field strength based on published simulation data, in a preliminary effort to determine the reliability of field strengths inferred to date in studies of galaxy cluster formation.

While this is not an easy task without actual simulation data on hand to analyze, the analysis work described in Miniati et al. (2000), based partly on ΛCDM simulations described in Cen & Ostriker (1994), provides enough information for us to make a rough estimate of the typical field strength. The simulation in question has the properties ${{H}_{0}}=60$ km s⁻¹ Mpc⁻¹, ${{{\Omega }}_{M}}=0.45$, ${{{\Omega }}_{{\Lambda }}}=0.55$, ${{{\Omega }}_{b}}=0.043$. We now give a relatively detailed description of our analysis of the information in Miniati et al. (2000), in order to make clear both the basis for the estimated field strength and the considerable uncertainty that attend these estimates.

Our starting point is Equation (33), which gives the postshock field strength due to the Biermann effect, assuming an unmagnetized upstream fluid. In order to use this equation to estimate field strengths, we need typical values of the shock speed D, the compression ratio ${{\rho }_{d}}/{{\rho }_{u}}$, and the tangential gradient $|{{\nabla }_{\bot }}({{k}_{B}}{{T}_{e}})|\equiv |\nabla ({{k}_{B}}{{T}_{e}})\times {\boldsymbol{n}} |$ that arise at primordial epochs.

For the sake of exploiting the information available in Miniati et al. (2000), we will choose z = 3 as our fiducial primordial epoch. Examining the top right panel of Figure 8 in Miniati et al. (2000), we conclude that a not-atypical value of $|{{\nabla }_{\bot }}({{k}_{B}}{{T}_{e}})|$ is

$$\begin{eqnarray}&&|{{\nabla }_{\bot }}({{k}_{B}}{{T}_{e}})|\approx \frac{{{k}_{B}}\times {{10}^{6}}\;{\rm K}}{5\;{\rm Mpc}}\approx {{10}^{-35}}\;{\rm erg}\;{\rm c}{{{\rm m}}^{-1}}.\end{eqnarray} \tag{ 57 }$$

Neglecting momentarily the fact that the simulation certainly fails to resolve the thermal precursor length scale ${{\lambda }_{T}}$ (Equation (60) below), the uncertainty in this estimate seems to be a factor of 2 or so.

Miniati et al. (2000) supply shock speeds at redshift z = 0 (see their Figure 5). These may be approximately shifted to z = 3 using Figure 10 of Miniati et al. (2000), which shows that the kinetic energy processed by shocks increased by a factor of 15 between z = 3 and z = 0. This suggests that temperatures increased by about that much (neglecting shock filling factor differences between the two redshifts) and that velocities increased by about a factor of 4. We had typical shock temperatures of $\sim {{10}^{6}}$ K at z = 3 in the data leading to the estimate in Equation (57). This corresponds to a temperature $\sim {{10}^{7}}$ K at z = 0. From the top right panel of Figure 5 of Miniati et al. (2000), this corresponds to a shock speed of about $6\times {{10}^{7}}$ cm s⁻¹ at z = 0, and hence $D\sim {{10}^{7}}$ cm s⁻¹ at z = 3.

Putting these numbers in Equation (33) and assuming the strong-shock limit ${{\rho }_{d}}/{{\rho }_{u}}=4$ appropriate to a $\gamma =5/3$ ideal gas, we obtain

$$\begin{eqnarray}\begin{array}{rcl} B & = & 2.9\times {{10}^{-22}}{\rm G}\times {{\left( \frac{D}{{{10}^{7}}\;{\rm cm\ }{{{\rm s}}^{-1}}} \right)}^{-1}} \\ {} & {} & \times \left( \frac{|{{\nabla }_{\bot }}{{k}_{B}}{{T}_{e}}|}{{{10}^{-35}}\;{\rm erg\ c}{{{\rm m}}^{-1}}} \right). \\ \end{array}\end{eqnarray} \tag{ 58 }$$

This value should be regarded as a lower limit since, as remarked above, the simulation does not resolve ${{\lambda }_{T}}$. An upper bound on B can be obtained by replacing, in Equation (57), the 5 Mpc length scale estimated from Figure 8 of Miniati et al. (2000) with ${{\lambda }_{T}}$ from Equation (60). This gives $B\lesssim {{10}^{-19}}$ G, i.e., a value that is a factor of about 10³ higher. Obviously, this is a very conservative bound, as it assumes temperature fluctuations of order unity over the diffusive scale ${{\lambda }_{T}}$.

It is clear from the very uncertain nature of the factors used to construct this estimate that the field strength given in Equation (58) is itself subject to considerable uncertainty. A few actual MHD simulations with a correct implementation of the Biermann effect seem required to establish how much field strength is in fact made available by the Biermann effect, for turbulent dynamo effects to amplify. Such simulations would be challenging given the spatial scale ${{\lambda }_{T}}$ that needs to be resolved.

5.1.1. Galaxy Cluster Formation Length Scales

We now estimate the physical length scales of the resistive magnetic precursor and the electron thermal precursor for the case of galaxy cluster formation using the physical parameters inferred in Section 5.1 from Miniati et al. (2000) for galaxy cluster formation at $z\simeq 3$. As in Section 5.1, we use ${{T}_{e}}={{10}^{6}}$ K as our reference temperature and $D={{10}^{7}}$ cm s⁻¹ as our reference shock speed. The assumed cosmological parameter ${{{\Omega }}_{b}}=0.043$ corresponds to a proton number density ${{n}_{i}}\approx {{10}^{-7}}$ cm⁻³ at z = 0, and hence to ${{n}_{i}}\approx 6\times {{10}^{-6}}$ at z = 3, which we take as the density upstream of the shocks. With these parameters, the Coulomb logarithm is ${\rm ln} {\Lambda }\approx 40$. Setting Z = 1, A = 1, we then have

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{B}} & = & 8.20\times {{10}^{-3}}{\rm cm}\times \frac{{\rm ln} {\Lambda }}{40} \\ {} & {} & \times \left( \frac{D}{{{10}^{7}}\;{\rm cm}\;{{{\rm s}}^{-1}}} \right)\times {{\left( \frac{{{T}_{e}}}{{{10}^{6}}\;{\rm K}} \right)}^{-3/2}}, \\ \end{array}\end{eqnarray} \tag{ 59 }$$

and

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{T}} & = & 5.01\;{\rm kpc}\times \xi \times \frac{{\rm ln} {\Lambda }}{40}\times \left( \frac{D}{{{10}^{7}}\;{\rm cm}{{{\rm s}}^{-1}}} \right) \\ {} & {} & \times {{\left( \frac{{{T}_{e}}}{{{10}^{6}}\;{\rm K}} \right)}^{5/2}}\times {{\left( \frac{{{n}_{i}}}{6\times {{10}^{-6}}{\rm c}{{{\rm m}}^{-3}}} \right)}^{-1}}. \\ \end{array}\end{eqnarray} \tag{ 60 }$$

The much more tenuous and considerably hotter primordial plasma reverses the relative sizes of ${{\lambda }_{B}}$ and ${{\lambda }_{T}}$ compared to the HEDP case considered above: ${{\lambda }_{B}}$ is now negligible, while ${{\lambda }_{T}}$ is an astrophysically significant length. It is worth comparing ${{\lambda }_{T}}$ to ${{\lambda }_{e}}$ and ${{\lambda }_{i}}$, the expected Spitzer mean free paths of electrons and ions in the plasma, as a consistency check. This is the product of the mean collision time τ with the thermal speed ${{(3{{k}_{B}}T/m)}^{1/2}}$. Using the expression for ${{\tau }_{e}}$ from Huba (2007), we have

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{e}} & = & \frac{3m_{e}^{1/2}{{\left( {{k}_{B}}{{T}_{e}} \right)}^{3/2}}}{4\sqrt{2\pi }{{n}_{e}}{\rm ln} {\Lambda }{{e}^{4}}}\times {{\left( \frac{3{{k}_{B}}{{T}_{e}}}{{{m}_{e}}} \right)}^{1/2}} \\ {} & = & \frac{3\sqrt{3}}{4\sqrt{2\pi }}\frac{{{\left( {{k}_{B}}{{T}_{e}} \right)}^{2}}}{{{n}_{e}}{\rm ln} {\Lambda }{{e}^{4}}} \\ {} & = & 2.50\times {{10}^{-1}}\;{\rm kpc}{{\left( \frac{{{T}_{e}}}{{{10}^{6}}\;{\rm K}} \right)}^{2}} \\ {} & {} & \times {{\left( \frac{{\rm ln} {\Lambda }}{40} \right)}^{-1}}\times {{\left( \frac{{{n}_{e}}}{6\times {{10}^{-6}}\;{\rm c}{{{\rm m}}^{-3}}} \right)}^{-1}}. \\ \end{array}\end{eqnarray} \tag{ 61 }$$

The mean free path is independent of particle mass, and for Z = 1 it is the case that ${{\lambda }_{e}}={{\lambda }_{i}}$. This length is the characteristic size of the viscous shock sheath and is reassuringly shorter than the semihydrodynamically scaled ${{\lambda }_{T}}$.

It is also worth considering whether the plasma is collisional, as we have implicitly assumed by using Spitzer-type transport coefficients. The magnetic field strength estimated in Section 5.1 gives rise to electron gyroradii ${{\lambda }_{c}}$ of characteristic size

$$\begin{eqnarray}\begin{array}{rcl} {{\lambda }_{c}} & = & {{\left( \frac{3{{k}_{B}}{{T}_{e}}}{{{m}_{e}}} \right)}^{1/2}}{{\left( \frac{eB}{mc} \right)}^{-1}} \\ {} & = & 41.4\;{\rm kpc}{{\left( \frac{{{T}_{e}}}{{{10}^{6}}\;{\rm K}} \right)}^{1/2}}\times {{\left( \frac{B}{3\times {{10}^{-22}}\;{\rm G}} \right)}^{-1}}. \\ \end{array}\end{eqnarray} \tag{ 62 }$$

Since ${{\lambda }_{c}}\sim 100{{\lambda }_{e}}$, the flow is comfortably collisional for the chosen physical parameters. An increase in the estimate of B by two orders of magnitude or more would change this conclusion; however, it should be noted that at least at the outset of the process of cosmic magnetogenesis envisioned here, field strengths were probably small enough that the collisional plasma approach is valid. Primordial field strengths that may have originated in early-universe phase transitions or in inflationary scenarios are poorly constrained by theory and observation, but they are believed to have plausible values $B\lt {{10}^{-22}}$ G in all but a few scenarios (see Widrow 2002; Widrow et al. 2012). It follows that, initially at least, the Spitzer-type transport coefficients that we have employed here are likely to be appropriate. This is in contrast to the apparent suppression of electron conductivity at later times (see, e.g., Markevitch et al. 2003; Russell et al. 2012; Sanders et al. 2013; ZuHone et al. 2013).

We have implemented the corrected Biermann effect algorithm, described in Section 3.4, within the publically available FLASH hydrodynamic simulation framework (Fryxell et al. 2000; Dubey et al. 2009; Tzeferacos et al. 2014). FLASH is a modular and extensible multiphysics scientific simulation software package that has been widely used for reactive compressible flows typical of astrophysical situations, HEDP applications, cosmology, computational fluid dynamics, and fluid–structure interactions. In particular, FLASH makes available both a two-temperature single-fluid model and resistive MHD, which makes it an ideal platform for testing the proposed algorithm.

Below we furnish an outline of the workings of the code, including the newly implemented Biermann effect algorithm. Much more detail about FLASH is available in the FLASH Users Guide.¹ We plan to include an implementation of the new Biermann effect algorithm in a future release of the code.

The two-temperature MHD model that we use for our numerical tests is expressed in conservation form by the following dynamical system:

$$\begin{eqnarray}&&\frac{\partial \rho }{\partial t}+\nabla \cdot \left( \rho {\boldsymbol{u}} \right)=0\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&\frac{\partial \rho {\boldsymbol{u}} }{\partial t}+\nabla \cdot \left[ \rho {\boldsymbol{uu}} -\frac{1}{4\pi }{\boldsymbol{BB}} +{\boldsymbol{1}} \left( P+\frac{{{{\boldsymbol{B}} }^{2}}}{8\pi } \right) \right]=0\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial }{\partial t}\left[ \rho \left( \frac{1}{2}{{{\boldsymbol{u}} }^{2}}+{{\epsilon }_{T}} \right)+\frac{{{{\boldsymbol{B}} }^{2}}}{8\pi } \right]+\nabla \cdot \left[ \rho {\boldsymbol{u}} \left( \frac{1}{2}{{{\boldsymbol{u}} }^{2}}+{{\epsilon }_{T}}+\frac{P}{\rho } \right) \right. \\ \quad +\frac{1}{4\pi }\left( -{\boldsymbol{u}} \times {\boldsymbol{B}} +\frac{{{c}^{2}}\eta }{4\pi }\nabla \times {\boldsymbol{B}} \right)\times {\boldsymbol{B}} \\ \quad +\left. \frac{c{{k}_{B}}}{4\pi e}{\rm ln} {{P}_{e}}\left[ \nabla \times \left( {{T}_{e}}{\boldsymbol{B}} \right) \right] \right]=0 \\ \end{array}\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&\begin{array}{l} \frac{\partial \rho {{s}_{e}}}{\partial t}+\nabla \left( \rho {\boldsymbol{u}} {{s}_{e}} \right)=-T_{e}^{-1}\nabla \cdot \left( -{{\kappa }_{e}}\nabla {{T}_{e}} \right) \\ \quad +\frac{\rho {{c}_{v,e}}}{{{T}_{e}}{{\tau }_{ei}}}({{T}_{i}}-{{T}_{e}}) \\ \end{array}\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}+\nabla \cdot \left[ {\boldsymbol{uB}} -{\boldsymbol{Bu}} -\frac{{{c}^{2}}\eta }{4\pi }\nabla {\boldsymbol{B}} +{{{\boldsymbol{G}} }^{(B)}} \right]=0.\end{eqnarray} \tag{ 67 }$$

These equations differ from Equations (8)–(12) only in that the Biermann term in Equation (65) replaces the one in Equation (10) to reflect Equation (46), while the Biermann term in Equation (67) replaces the one in Equation (12) to reflect Equation (40).

FLASH advances the MHD equations using an unsplit staggered mesh (USM) scheme that prevents the development of magnetic monopoles due to numerical noise in the induction equation, and which is described in Lee (2013).

The thermal conduction and heat-exchange source terms on the right-hand side of Equation (66) are computed in a time-split manner, separately from the ideal MHD advance. In particular, thermal conduction advance is fully implicit and works as described in Section 17.1.4 of the FLASH Users Guide, while heat exchange operates as described in Section 16.5 of the FLASH Users Guide. The remaining, advective part of Equation (66) is implemented by treating s_e as a passively advected mass scalar. This equation requires an EOS implementation capable of using electron entropy as an input variable. For the current case of perfect gas EOS and total ionization, the standard Sackur–Tetrode equation for entropy is implemented in the EOS.

The resistive terms in Equations (65) and (67) are not treated implicitly, but rather are added directly to the MHD fluxes and are thus treated explicitly, placing diffusive stability limitations on the maximum time step. The Biermann terms in Equations (65) and (67) are also added explicitly to the MHD fluxes. These flux modifications are performed in such a way as to preserve the solenoidal character of the magnetic field, by adding them to the ideal MHD Godunov fluxes before these are interpolated to the edge-centered electric fields of the USM scheme, after which each face-centered normal magnetic field component is updated by the circulation integral of the electric field along the edges bounding the face (Lee 2013).

We have not performed an analysis of the time step limitation imposed by the Biermann effect, trusting instead that the small magnitude of the effect in the cases considered makes it unlikely that Biermann time step constraints could dominate those due to the hyperbolic and diffusive terms. In this connection, it is noteworthy that the Biermann term, being quadratic in fluid derivative terms, makes no contribution to wave dispersion relations obtained by linearization about uniform solutions and therefore does not appear to affect ordinary wave speeds. Intuitively, one would therefore expect any time step limitations due to the Biermann effect to be higher-order corrections, hardly affecting the numerical evolution of plasmas for which the Biermann term is small in magnitude. Nonetheless, a full analysis of the effect of the Biermann term on the hyperbolic structure of the system of PDEs would be useful—possibly even necessary—for cases of intense field generation. Such an analysis would be somewhat complicated by the non-quasi-linear structure of the Biermann term, necessitating an extension of the PDE system to quasi-linear normal form.

The simulations presented below represent idealized situations simplified for the sake of verifying the code and for illustrating the numerical and physical principles under discussion. In particular, the conductivity ${{\kappa }_{e}}$, electron–ion equilibration timescale ${{\tau }_{ei}}$, and resistivity η are all treated as constants. FLASH does have the capability to use Spitzer-type functions of the thermal state for these parameters, but this would complicate the results unnecessarily without adding any real value to these verification tests.

In what follows, we assume fully ionized hydrogen—A = 1, Z = 1, adiabatic index $\gamma =5/3$, ${{c}_{v,e}}=\frac{3}{2}{{k}_{B}}{{N}_{A}}$, where N_A is Avogadro's number and k_B is Boltzmann's constant. Simulations are conducted in cylindrical coordinates, assuming azimuthal symmetry, and are therefore two-dimensional. In addition, we impose a reflection boundary on the R-axis (where R is the perpendicular distance from the axis of rotational symmetry, which is to say, R is the cylindrical radius) so that the domain represents a hemisphere of a solution with reflection symmetry about that axis. In every case, the domain is a 4 cm radius cylinder that extends 4 cm in the z direction from the R axis. The boundary conditions at R = 4 cm and at z = 4 cm are outflow.

All of the verification tests that we present are variations on the theme of a Sedov-esque explosion: an analytic self-similar Sedov solution is (1) modified to prevent temperatures from diverging at the center by setting all variables constant inside some chosen radius; (2) smoothed with a Gaussian near the shock, to ease instabilities that can otherwise appear near the shock; (3) where required, distorted to an ellipsoidal profile capable of producing real Biermann effect fields; and (4) scaled to velocities and pressures estimated to produce usable simulation times given the resolutions and domain size in play. Values of shock ellipticity, ${{\kappa }_{e}},$and η are then chosen to produce reasonable thermal and magnetic precursor region sizes and reasonable field generation rates, and values of ${{\tau }_{ei}}$ are chosen to be harmless. In every case, we run an initial model for a time without the Biermann effect, so as to allow transients in the two-temperature hydrodynamic variables to subside. Then we restart the calculation with Biermann field generation turned on.

Magnetic fields are always azimuthal in these verification tests—since the gradients of T_e and P_e are always in the $R-z$ plane, the Biermann effect only generates a nonzero field along the azimuthal direction.

6.2.1. Null Test: A Spherical Shock

It is not a simple matter to construct a nontrivial analytic verification solution of a plasma generating magnetic field through the Biermann effect. A trivial solution, on the other hand, may be straightforwardly constructed by imposing symmetry requirements that align the gradients $\nabla {{P}_{e}}$ and $\nabla {{T}_{e}}$, thus guaranteeing zero-field generation. We use a spherical Sedov-like explosion as an example of such a test. This is in fact the test that revealed the Biermann catastrophe in the first place (Fatenejad et al. 2013).

We assume an initial shock radius of 1.15 cm, in an ambient medium with $\rho =1$ g cm⁻³ and P set to a negligible value. We choose an initial velocity scale inside the shock that leads to a shock velocity that is on average about $3.2\times {{10}^{4}}$ cm s⁻¹ during the course of the simulation. Other physical parameters are ${{\kappa }_{e}}=1.0\times {{10}^{12}}$ erg K ⁻¹ cm⁻¹ s⁻¹ and ${{\tau }_{ei}}=2\times {{10}^{-14}}$ s. The initial solution is advanced for 10⁻⁵ s with no field generation and then restarted and run for another 10⁻⁵ s with field generation turned on. We repeated this experiment with the correct Biermann flux term of Equation (40), with the naive flux of Equation (41), and with the time-split source term of Equation (15), so as to compare the performance of the three algorithms.

The results of this study are displayed in Figure 2. The top left panel is a color map view of the electron temperature T_e. The dashed shock-crossing line across the bottom of this figure is the location of the data displayed in the top right panel, which shows electron and ion temperatures. The shock is visible as the peak in T_i. The thermal precursor is clearly visible, as is the continuous behavior of T_e, which changes slope at the shock in accordance with standard plasma shock theory (Zel'dovich & Raizer 2002, Chapter 7, Section 12). The ion temperature also rises in the precursor region in consequence of the electron–ion equilibration term.

Zoom In Zoom Out Reset image size

The bottom left panel of Figure 2 shows the magnetic field strength generated by the correct treatment of Equation (40). There is a spurious nonzero field due to numerical noise that is clearly being generated. However, the bottom right panel shows that this field is converging to zero with increasing resolution. This figure displays the square root of the total magnetic energy in the domain. Since the analytic solution for the magnetic field strength is zero everywhere, this quantity functions as an un-normalized L²-norm of the difference between the analytic and numerical solutions. The blue circles show the convergence with resolution of the correct formulation. This convergence behavior is in contrast to the behavior of the other two algorithms, which fail altogether to converge.

6.2.2. Verification with an Ellipsoidal Shock

Next we exhibit simulations designed to produce nonzero values of the magnetic field strength. We start with an ellipsoidal shock surface, obtained by distorting the spherical Sedov solution. The initial configuration has a semimajor axis (aligned with the R-axis) of 2 cm and a semiminor axis (aligned with the z-axis) of 1.667 cm. The ambient medium is uniform with $\rho =1$ g cm⁻³ and P set to a negligible value. We choose a velocity scale that leads to shock velocities of about $6\times {{10}^{6}}$ cm s⁻¹. The higher shock speed is chosen to produce somewhat intense magnetic field strengths. To keep the thermal precursor resolved, we increase the conductivity to ${{\kappa }_{e}}=1.0\times {{10}^{13}}$ erg K⁻¹ cm ⁻¹ s⁻¹. We also set ${{\tau }_{ei}}={{10}^{-14}}$ s. The initial solution is advanced for $4.86\times {{10}^{-8}}$ s with no field generation and then restarted and advanced to a simulation time $t=3\times {{10}^{-7}}$ s with field generation turned on.

The advance of the shock during the period when magnetic field is generated is illustrated by the two pressure color map plots in the top panels of Figure 3. The bottom left panel of Figure 3 shows the magnetic field intensity generated according to the correct flux of Equation (40), ranging into the tens of gauss. The solid red line in the figure shows the shock location, while the red arrows display the shock unit normal vector. The bottom right panel displays the difference between the correct and the naive flux formulations, which is substantial at the outgoing shock surface.

Zoom In Zoom Out Reset image size

We do not have an analytical solution for the magnetic field distribution to compare to the output of these simulations. We do, however, have the relation between shock quantities expressed by Equation (32), which determines the jump condition for ${\boldsymbol{B}}$. By locating points adjoining the shock and computing the local shock velocity at those points, we can then verify that Equation (32) in fact obtains, to some accuracy limited by the numerical approximation. This is, in principle, a nontrivial verification, since the code does not know about the jump condition Equation (32) as such—it only knows the fluxes expressed by Equation (40), which were inferred from the jump condition. Recovery of the jump condition thus constitutes a verification that the code correctly models the theory.

To perform this verification, we first locate cells adjoining the shock by the method described in the Appendix, which yields a one-cell-wide shock surface by fitting the mass, momentum, and energy R–H conditions to the neighboring data, using a speed-of-sound weighted inner product on the state space, and treating the shock speed D as a free fit parameter. The shocked cells are the ones that fit the R–H conditions with small fit residuals (${{R}^{2}}\lt 0.01$), together with some other natural auxiliary consistency conditions described in the Appendix. The fitted shock speed is then used in the verification of Equation (32).

In the absence of normal component field B_n, Equation (32) becomes

$$\begin{eqnarray}\begin{array}{lll} {} & {} & (D-{{u}_{d}}){{B}_{d}}-(D-{{u}_{u}}){{B}_{u}} \\ {} & {} & +\frac{c{{k}_{B}}}{e}\left( {{n}_{z}}\frac{\partial {{T}_{e}}}{\partial R}-{{n}_{R}}\frac{\partial {{T}_{e}}}{\partial z} \right)\left[ {\rm ln} {{P}_{e,d}}-{\rm ln} {{P}_{e,u}} \right]\;\equiv \\ {} & {} & {{a}_{d}}-{{a}_{u}}+{{b}_{d}}-{{b}_{u}}=0, \\ \end{array}\end{eqnarray} \tag{ 68 }$$

where we have defined "advection" terms ${{a}_{u,d}}$ and "Biermann" terms ${{b}_{u,d}}$ by

$$\begin{eqnarray}&&{{a}_{d}}\equiv (D-{{u}_{d}}){{B}_{d}}\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&{{a}_{u}}\equiv (D-{{u}_{u}}){{B}_{u}}\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{{b}_{d}}\equiv \frac{c{{k}_{B}}}{e}\left( {{n}_{z}}\frac{\partial {{T}_{e}}}{\partial R}-{{n}_{R}}\frac{\partial {{T}_{e}}}{\partial z} \right){\rm ln} {{P}_{d}}\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&{{b}_{u}}\equiv \frac{c{{k}_{B}}}{e}\left( {{n}_{z}}\frac{\partial {{T}_{e}}}{\partial R}-{{n}_{R}}\frac{\partial {{T}_{e}}}{\partial z} \right){\rm ln} {{P}_{u}}.\end{eqnarray} \tag{ 72 }$$

The sum of terms in Equation (68) is required to be zero. In a discretized numerical PDE integration, this really means that the terms must cancel up to some truncation or rounding precision, which is expressed relative to the largest of the magnitudes in play. We therefore define the "shock condition parameter" C by

$$\begin{eqnarray}&&C\equiv \frac{{{a}_{d}}-{{a}_{u}}+{{b}_{d}}-{{b}_{u}}}{{\rm max} (|{{a}_{d}}|,|{{a}_{u}}|,|{{b}_{d}}|,|{{b}_{u}}|)}.\end{eqnarray} \tag{ 73 }$$

We calculate the value of C at cells along the shock front, at each of our four resolution levels, for both the correct flux and the "naive" flux implementations of the Biermann effect. We display cumulative distributions of C in the top panels of Figure 4. It is evident from these figures that the distribution of C is in fact centered near zero for both flux implementations. The widths of the distributions behave very differently, however. In the case of the correct Biermann flux implementation, the distributions get narrower with each refinement of resolution, providing some evidence of convergence with resolution to the expected result. In the case of the naive flux implementation, there is no such evidence of convergence. The same observations can be made of the plots in the middle panels of Figure 4, which summarize the means and sample standard deviations of the C-distributions, as a function of resolution. The convergence of the correct flux implementation and the convergence failure of the naive implementation are clear in these figures.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 4 show the evolution of total magnetic energy in the domain as a function of time, for the different resolutions and for the two flux implementations. The correct flux implementation appears to be converging (although not in any strong sense converged) by this measure, even at late time, whereas any evidence of convergence in magnetic field energy is simply lacking for the naive flux implementation.

6.2.3. The Resistive Magnetic Precursor

In this set of simulations, we maintain the simulation parameters described in Section 6.2.2 but also turn on the resistivity η to a finite positive value, $\eta =7.8\times {{10}^{5}}\;{{{\rm s}}^{-1}}$, chosen in conjunction with the other parameters to yield an easily discernible exponentially decaying resistive precursor described by Equation (52). We repeat the simulation strategy of Section 6.2.2, letting the simulation advance to $t=4.86\times {{10}^{-8}}$ s with no field generation, and then restarting it and advancing to a simulation time $t=3\times {{10}^{-7}}$ s with field generation (using the correct flux formulation) turned on. This advance time is adequate for the establishment of the expected precursor structure, as we demonstrate below.

The left panel of Figure 5 shows a color map of the distribution of magnetic field strength across the domain. The shock location is shown by the solid red line, and the superposed vectors indicate the shock-normal direction. The magnetic field is evidently smoothed out by the resistivity, as can be seen by a direct comparison with the bottom left panel of Figure 3. The precursor is already evident in this figure, as the substantial amount of magnetic field that has "leaked" ahead of the shock.

The diagonal white dashed shock-crossing line in the left panel of Figure 5 illustrates the line along which data were extracted to produce the right panel of Figure 5. This figure shows the magnetic field values plotted along that line. The vertical solid red line at X = 0 marks the location of the shock. The exponential decay of the field is easily visible in this figure.

At each shock location, the value of η can be combined with the local shock velocity to calculate ${{\lambda }_{{\rm Predicted}}}$, the expected decay length of the precursor, according to Equation (53). At the same time, the actual decay length ${{\lambda }_{{\rm Measured}}}$ can be measured directly by comparing the field strength at two suitably selected distances along the local shock normal. We have performed both calculations at each shock location and compared them. The results are displayed in Figure 6, where we plot the cumulative distribution of ${{\lambda }_{{\rm Measured}}}/{{\lambda }_{{\rm Predicted}}}$. The distribution appears to be consistent with a mean value of 1, as expected, with some scatter. The scatter is not surprising, since the expression in Equation (53) for ${{\lambda }_{{\rm Predicted}}}$ is an approximation, assuming, as it does, propagation of magnetic field from a planar shock. Since the shock is, of course, not planar, the value of the Biermann generation rate changes appreciably across the shock over distances not hugely different from ${{\lambda }_{{\rm Predicted}}}$ itself. Under the circumstances, then, the level of agreement between observation and prediction is satisfactory.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2.4. The Thermal Magnetic Precursor

In our final simulation, we explore the properties of the thermal magnetic precursor by returning to a nonresistive simulation, similar to the ones discussed in Section 6.2.2, and differing from them only in that ${{\kappa }_{e}}$ is 10 times larger: ${{\kappa }_{e}}=1.0\times {{10}^{14}}$ erg K⁻¹ cm⁻¹ s⁻¹. This broadens the thermal precursor zone to about 0.2 cm, making it easier to discern.

The results at the final time step are illustrated in Figure 7. The top left panel shows the distribution of magnetic field in the domain. The solid red line in the figure shows the location of the shock, while the arrows show the direction of the shock normal. In the top right panel we see the electron and ion temperatures along the shock-crossing line shown in the previous panel. The shock location coincides with the sharp drop in T_i, while the thermal conduction precursor zone, which has an exponential structure for the constant conductivity model used here, may be clearly seen in the upstream behavior of T_e.

Zoom In Zoom Out Reset image size

The magnetic field structure along the shock-crossing line is shown in the bottom left panel of the figure. The precursor is difficult to see in this plot, so the bottom right panel shows $|B|$ on a semi-log plot. The magnetic precursor structure is seen in this figure. It clearly has a much lower field intensity than the fields generated by the Biermann effect at the shock. We emphasize again, however, that in part this is due to the highly simplified constant-conductivity model adopted for this work. A true Spitzer-type conductivity, with a $T_{e}^{5/2}$ dependence, would create much sharper gradients at the upstream termination of the precursor zone, potentially generating more intense fields—relative to those generated at the shock—than those shown here.

The first condition above is tantamount to insisting that the shocked cell be sited at a position where the gradient of P is steepest. In general, volume-based hydrodynamic advance schemes spread shocks over several cells. If we were not to impose this condition, we would obtain a thickened shell of shocked cells, which would complicate the criterion (specified below) for classifying neighboring cells as adjoining the shock upstream or downstream.

The condition is very easy to enforce: in the immediate neighborhood of the cell being tested for "shockedness," we check the adjoining cells closest to the direction of the normal vector ${\boldsymbol{n}}$ and its mirror image $-{\boldsymbol{n}}$. This is done by forming the vector of integers ${\Delta }{\boldsymbol{i}} \equiv {\rm NINT}[{\boldsymbol{n}} /{\rm max} ({{n}_{1}},{{n}_{2}},{{n}_{3}})]$ (where the NINT function ascribes the nearest integer to each component of its real vector argument) and adding it to the cell discrete index vector ${\boldsymbol{i}}$, to reach the cells at ${\boldsymbol{i}} \pm {\Delta }{\boldsymbol{i}}$. If a stencil-based estimate of $|\nabla P|$ is greater in the candidate cell than in the two so-chosen adjoining cells, then the condition is passed successfully.

An easy sanity check for "shockedness" of a cell is that cells downstream should unfailingly be (a) compressed, (b) accelerated (in the $+{\boldsymbol{n}}$-direction), and (c) heated, with respect to cells upstream.

We introduce a shock-width parameter h, such that the shock-adjoining cells "upstream" and "downstream" of the shocked cell are respectively located at $\pm h\times {\Delta }\times {\boldsymbol{n}}$ relative to the cell under study, where Δ is the length of a cell side. For FLASH with the HLL Riemann solver (Toro 2009) and second-order reconstruction, an appropriate value of h is 1.5. We then simply verify that the fluid variables ρ, ${\boldsymbol{u}}$, and e_T (the total thermal energy) satisfy the conditions $\left[ \rho \right]_{u}^{d}\gt 0$, $\left[ {\boldsymbol{u}} \cdot {\boldsymbol{n}} \right]_{u}^{d}\gt 0$, and $\left[ {{e}_{T}} \right]_{u}^{d}\gt 0$.

The R–H conditions, which express flux conservation in the frame of the shock, have the form

$$\begin{eqnarray}&&D\times ({{{\bf \Phi }}_{{\rm d}}}-{{{\bf \Phi }}_{{\rm u}}})={\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm d}}})-{\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm u}}}).\end{eqnarray} \tag{ A1 }$$

Here ${\bf \Phi }$ is the state vector of conserved field quantities, ${{{\bf \Phi }}^{T}}\equiv [\rho ,\rho {\boldsymbol{u}} \cdot {\boldsymbol{n}} ,\rho \mathcal{E}]$, that is, the mass, normal momentum, and energy densities. ${\boldsymbol{F}}$ denotes the vector of fluxes corresponding to ${\bf \Phi }$. The subscripts "d" and "u" denote "downstream" and "upstream" states, respectively (Toro 2009).

We fit R–H conditions to the upstream/downstream data defined using the shock-width parameter h above. In this fit, the shock speed D is a free parameter to be adjusted to minimize a fit residual. Given a positive-definite inner product $\left[ {{{\Phi }}_{1}},{{{\Phi }}_{2}} \right]$ on the vector state space in Equation (A1), we may define the normalized residual R² as

$$\begin{eqnarray}&&{{R}^{2}}\equiv \frac{[\left( D{\Delta }{\bf \Phi }-{\Delta }{\boldsymbol{F}} \right),\left( D{\Delta }{\bf \Phi }-{\Delta }{\boldsymbol{F}} \right)]}{[{\Delta }{\boldsymbol{F}} ,{\Delta }{\boldsymbol{F}} ]},\end{eqnarray} \tag{ A2 }$$

where ${\Delta }{\bf \Phi }\equiv {{{\bf \Phi }}_{d}}-{{{\bf \Phi }}_{u}}$ and ${\Delta }{\boldsymbol{F}} \equiv {\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm d}}})-{\boldsymbol{F}} ({{{\bf \Phi }}_{{\rm u}}})$. We may then easily minimize R² with respect to D, obtaining

$$\begin{eqnarray}&&{{D}_{{\rm Min}}}=\frac{[{\Delta }{\bf \Phi },{\Delta }{\boldsymbol{F}} ]}{[{\Delta }{\bf \Phi },{\Delta }{\bf \Phi }]}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&R_{{\rm Min}}^{2}=1-\frac{{{[{\Delta }{\boldsymbol{F}} ,{\Delta }{\bf \Phi }]}^{2}}}{[{\Delta }{\boldsymbol{F}} ,{\Delta }{\boldsymbol{F}} ][{\Delta }{\bf \Phi },{\Delta }{\bf \Phi }]}.\end{eqnarray} \tag{ A4 }$$

Given these definitions, we consider that the R–H conditions are satisfied if R²_Min is less than some small threshold value.

It remains to define the inner product $[\cdot ,\cdot ]$ used in these expressions. It is obvious that the naive unweighted sum-of-squares of the components of vectors such as ${\Delta }{\bf \Phi }$ is not an acceptable inner product, as it is dimensionally senseless. We require at a minimum a positive-definite metric that brings all vector components to the same physical dimensions, so that we do not wind up adding mass densities to energy densities, and so on. An additional desirable requirement is that all three components of Equation (A1) should contribute similar magnitudes to Equations (A3) and (A4), so that they are all weighted equally in the fit.

A simple metric that accomplishes these requirements may be built using the local sound speed c_s of the candidate shocked cell:

$$\begin{eqnarray}&&\left[ {{{\bf \Phi }}_{1}},{{{\bf \Phi }}_{2}} \right]\equiv {\bf \Phi }_{1}^{T}{\boldsymbol{M}} {{{\bf \Phi }}_{2}},\end{eqnarray} \tag{ A5 }$$

where

$$\begin{eqnarray}{\boldsymbol{M}} \equiv \left[ \begin{array}{ccccccccccccccc} 1 & 0 & 0 \\ 0 & c_{s}^{-2} & 0 \\ 0 & 0 & c_{s}^{-4} \\ \end{array} \right].\end{eqnarray} \tag{ A6 }$$

This choice produces compatible dimensions and comparable magnitudes, because in and downstream of the shock we have ${{c}_{s}}\sim {{({{k}_{B}}T/{{m}_{i}})}^{1/2}}\sim D$, so that in the frame of the shock, the mass flux is $\rho D\sim \rho {{c}_{s}}$, the momentum flux has magnitude $P\sim \rho c_{s}^{2}$, and the energy flux has magnitude $\rho {{D}^{3}}\sim \rho c_{s}^{3}$.

We adopt this metric in computing Equations (A3) and (A4). We find that a threshold of $R_{{\rm Min}}^{2}\lt 0.01$ is adequate for identifying cells satisfying the R–H conditions with few false positives and false negatives, in the simulations exhibited in this paper.

Our final criterion is easily stated and checked: shocks are supersonic upstream and subsonic downstream, so that characteristics from the family implicated in the shock converge on the shock (see Courant & Friedrichs 1948). This is an essential stability criterion for the solution. Using the shock speed D_Min calculated while fitting the R–H conditions to the data, this condition is simply

$$\begin{eqnarray}&&{{{\boldsymbol{u}} }_{d}}\cdot {\boldsymbol{n}} +{{c}_{s,d}}\gt {{D}_{{\rm Min}}}\gt {{{\boldsymbol{u}} }_{u}}\cdot {\boldsymbol{n}} +{{c}_{s,u}},\end{eqnarray} \tag{ A7 }$$

where ${{c}_{s,d}}$ and ${{c}_{s,u}}$ are the downstream and upstream sound speeds, respectively, as determined with respect to the shock width h defined above.